% These scripts were made to get
% scalar mag field anomaly data for
% VM Tiwari
%India
latmin = -89.5;
latmax = 89.5;
logmin = 0.5;
logmax = 359.5;
incr = 1;
% AFrica
% latmin = -30;
% latmax = 30;
% logmin = 0;
% logmax = 60;
% incr = 0.5;


tha = 90-(latmin:incr:latmax); % equal area 
phi = logmin:incr:logmax;
fday = 1800;
alt = 0;


% getting the main field (internal part) B[x,y,z,f]
pomme25_cf = 'C:\Manoj\ocean\pomme\pomme-2.5.cof';
sm2geo_cf  = 'C:\Manoj\ocean\pomme\sm2geo_coeff';
gsm2geo_cf = 'C:\Manoj\ocean\pomme\gsm2geo_coeff';

imf_by  = -5.0;
e_st    = -10.0;
i_st    = 0.3;

model  = pomme_init(pomme25_cf, sm2geo_cf, gsm2geo_cf);

pomme_control(model,'DEGREE RESET');
pomme_control(model,'DEGREE 16 90');
%     % -90 lat south pole so 90- (-90) = 180
%     % -180 long 
%     %% SINGLE CALCULATION, RETURNS:  bm_internal(X,Y,Z,F) and be_external(X,Y,Z)
%     flag = 1;
disp('Computing fields...............');

for i = 1:length(phi),
lat = tha;
on = ones([1,length(tha)]);
[c,a]=pomme_calc(model,fday*on,(6371.2+alt)*on,lat*pi/180,on*phi(i)*pi/180,on*e_st,on*i_st,on*imf_by);
B(i,:,:) = c';
end;

figure(2);
worldmap([latmin,latmax],[logmin,logmax]);
[c,ch] = contourfm(90-tha,phi,squeeze(B(:,:,3))',10);

for i = 1:4,
         B1(i,:,:) = B(:,:,i);
     end;
      clear B;
B=B1;
clear B1;


% Computing the main field to project the anomalies in main field direction
pomme_control(model,'DEGREE RESET');
pomme_control(model,'DEGREE 1 16');
%     % -90 lat south pole so 90- (-90) = 180
%     % -180 long 
%     %% SINGLE CALCULATION, RETURNS:  bm_internal(X,Y,Z,F) and be_external(X,Y,Z)
%     flag = 1;
disp('Computing fields...............');

for i = 1:length(phi),
lat = tha;
on = ones([1,length(tha)]);
[c,a]=pomme_calc(model,fday*on,(6371.2+alt)*on,lat*pi/180,on*phi(i)*pi/180,on*e_st,on*i_st,on*imf_by);
B2(i,:,:) = c';
end;
 
for i = 1:4,
         B1(i,:,:) = B2(:,:,i);
end;
      
bscn = dot(squeeze(B(1:3,:,:)),squeeze(B1(1:3,:,:)),1)./(B1(4,:,:));
% 
% 
% % subplot(411);
worldmap([latmin,latmax],[logmin,logmax]);
[c,ch] = contourfm(90-tha,phi,squeeze(B(3,:,:))',10);
colorbar;
title('B_z');

% subplot(312);
% worldmap([latmin,latmax],[logmin,logmax]);
% [c,ch] = contourfm(90-tha,phi,squeeze(B(4,:,:))',10);
% colorbar;
% title('|F|');
% 
% subplot(313);
worldmap([latmin,latmax],[logmin,logmax]);
[c,ch] = contourfm(repmat((90-tha)',[1,length(phi)]),repmat(phi,[length(tha),1]),squeeze(bscn)',10);
title('Scalar Anomaly');
colorbar;

lat = repmat((90-tha)',[1,length(phi)]);
long = repmat(phi,[length(tha),1]);
Sc_Pom_25 = squeeze(bscn)';
Bz_Pom_25 = squeeze(B(3,:,:))';
